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Summary. Using the concept developed in earlier papers, that 
the cosmic rays originate in three different main sites, a) the 
normal supernova explosions into the interstellar medium, b) 
the supernova explosions into a stellar wind, and c) powerful 
radio galaxies, we demonstrate in this paper that the spec- 
trum and chemical abundances above 10 4 GeV can be well 
understood. Using existing data on the chemical composition 
of cosmic rays near TeV energies as a constraint, we adjust 
the parameters of the model to fit the shower size data from 
the Akeno experiment; this is necessary since the derivation 
of an all particle spectrum involves an assumption about the 
chemical composition of the cosmic rays and so we have to 
fit the shower size data first and then derive the all particle 
spectrum. We present a successful fit to the shower size data 
which allows us to draw three main conclusions: a) For most 
of the energy range above 10 5 GeV the wind explosions can 
account for both chemical composition and spectrum includ- 
ing the knee feature, b) the highest particle energies required 
from the stellar wind explosions imply a magnetic field in the 
preexisting stellar wind of at least 3 Gauss at a fiducial dis- 
tance of 10 14 cm, c) the chemical abundances above 10 5 GeV 
are dominated by heavy nuclei such as Neon and higher. The 
bump observed in the all particle spectrum below the knee gets 
weakened with the proper treatment of the energy dependent 
chemical composition. At the high energy end we obtain an 
estimate of the extragalactic flux of protons. 



Key words: 

Shockwaves 



Cosmic Rays - Plasma Physics - Supernovae 



1. Introduction 

With the introduction of a new concept to treat particle ac- 
celeration in shocks where the shock normal is perpendicular 
to the prevailing magnetic field, Biermann (1993, paper CR I) 
has been able to interpret the cosmic ray spectrum from GeV 
to EeV energies, including the exact spectrum and the feature 
of the knee. The cosmic ray spectrum is composed of three 
components: 
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1) The explosions of normal supernovae into an approx- 
imately homogeneous interstellar medium drive blast waves 
which can accelerate to about 10 5 GeV for Hydrogen (La- 
gage and Cesarsky 1983). For these particles the spectrum is 
near —2.75 after taking leakage from the Galaxy into account. 
Particles get accelerated continuously during the expansion 
of the spherical Shockwave, with the highest particle energy 
reached at the beginning of the adiabatic expansion, the Sedov 
phase. However, after acceleration particles lose energy due to 
adiabatic expansion and hence, the maximum particle energy 
which is relevant, is the maximum particle energy achieved at 
the end of the Sedov phase, when cooling sets in and the shell 
can break up. At that point the accelerated particle population 
gets mixed with the interstellar medium. We will refer to this 
phase and the associated particle population as Sedov phase 
explosion and Sedov phase acceleration subsequently. We note 
that the strongest dependence of this maximum energy is on 
the density of the interstellar medium with the highest particle 
energy being achieved in the most tenuous medium. Protons 
can get accelerated up to about 10 5 GeV in the tenuous part 
of the interstellar medium. A detailed discussion of the accel- 
eration process and the differences between protons, heavier 
nuclei and electrons as well as the observed radio emission 
from normal supernova remnants is made in Biermann and 
Strom (1993, paper CR III). 

2) Explosions of stars into their former stellar winds (like 
Wolf Rayet stars) produce particles up to energies of about 
310 9 GeV (for iron, and rigidity dependent); these particles 
have a slightly flatter spectrum of near —2.67 (this is 0.08 flat- 
ter than the Sedov phase explosion spectrum) up to a rigidity 
dependent bend in the spectrum and then beyond that a 
steeper spectrum of near —2.97 up to the rigidity dependent 
cutoff. This spectral index difference is due to a difference in 
the detailed acceleration efficiency, which in turn depends on 
particle drift energy gains. Below the knee energy the drifts are 
stronger because the particles experience a stronger curvature 
drift due to turbulence, and beyond the knee energy the par- 
ticles just experience the basic curvature and gradient of the 
stellar wind, since their Larmor radius gets large (the details 
are explained in paper CR I). The discussion (Biermann and 
Cassinelli 1993, paper CR II) of the acceleration of particles in 
the strong shocks in radiosupernovae on the one hand, and the 
weaker shocks in the winds of OB and WR stars on the other 
hand yields the results: a) the magnetic field in the winds of 
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WR stars is of order 3 Gauss at a fiducial distance of 10 cm 
from the star, b) that thus particle energies up to about 3 10 9 
GeV in iron nuclei are possible from acceleration in strong 
supernova shocks in such winds, and c) that there is a critical 
velocity of order 10 9 cm/sec below which electron injection is 
strongly inhibited. 

3) The hot spots of Fanaroff Riley class II radio galaxies 
produce particles with even higher energies, up to near 10 11 
GeV; their spectrum is approximately E~ 2 up to the pileup 
just below the cutoff due to the interaction with the cosmolog- 
ical microwave background. This extragalactic component has 
been modelled by Rachen and Biermann (1992, 1993: paper 
UHE CR I): this latter model has been shown to be consistent 
with Fly's Eye data on the proton contribution below EeV 
energies (Gaisser et al. 1993, Rachen, Stanev, Biermann 1993, 
paper UHE CR II). 

The main fraction of the cosmic rays above 10 TeV is ac- 
celerated in shocks that traverse stellar winds. Stellar wind 
magnetic fields have asymptotically a Parker type topology, 
with Bff, ~ 1/r dominating over most of 4ir and a polar cap 
where the radial field B r ~ 1/r 2 dominates. In CR I Biermann 
developed a theory to treat perpendicular shocks in a spher- 
ical geometry for particle acceleration. This process produces 
E~' 7 ^ particle spectra in the limiting case of large shock ve- 
locities. In the geometrically small polar cap region, however, 
the geometry is that of a parallel shock, which in the limit of 
large shock velocities produces E~ 2 particle spectra. Although 
the total amount of energy carried by the polar cap component 
is small, we shall show that it has important contribution in 
the region of the knee. 

The Parker spiral magnetic field topology requires a region 
near the pole where the magnetic field is almost exactly radial; 
further out from the pole the magnetic field becomes helical to 
approach asymptotically the tight Archimedian spiral at the 
equator. There is a critical angle not far from the pole where 
the Larmor radii of the particles accelerated in the pole region 
overlap with the Larmor radii of the particles accelerated in the 
Archimedian spiral region; this critical angle separates the two 
regions where the acceleration is well approximated by a shock 
normal parallel to the magnetic field yielding an E~ 2 spectrum 
for strong shocks, and the region where the acceleration is best 
described by shock normal perpendicular to the magnetic field 
giving an E~~ '/ 3 spectrum. Since the pole region where the 
parallel acceleration dominates is only about 0.01 of 4ir the 
contribution of this region to the overall spectrum is very 
small except for the particle energies near 10 6 GeV, i.e. near 
the bend in the overall spectrum and the identical cutoff for the 
polar cap particles; the polar cap particles have a steep cutoff 
at the bend energy because at that energy the acceleration 
is limited by lack of space as discussed above. Fig. 1 shows 
a a schematic representation of the composite model we are 
discussing, broken down into its four components. 

Here we wish to test the overall model proposed by ask- 
ing whether it can successfully account for the spectrum and 
chemical composition at particle energies beyond 10 4 GeV: the 
special difficulty in this endeavour is the fact that we know the 
chemical abundances near TeV particle energies already, and 
extrapolating these spectra does not give the bump and knee 
of the well established overall particle spectrum. On the other 
hand, the particle energies themselves and the known chemical 
abundances can indeed be understood with the concept that 



Fig. 1. The generic spectrum of one nuclear species with all four 
cosmic ray components, here for Hydrogen: Component 1 is due to 
the Sedov phase explosions in the interstellar medium, component 2 
is due to stellar winds, with 3 the polar cap component, and compo- 
nent 4 is extragalactic. E\ is the cutoff energy for component 1, and 
i?2 is the bend energy for the wind component and also the cutoff 
energy for the polar cap component. 



explosions into a wind like that of a WR star are their origin 
(Volk and Biermann 1988, Silberberg et al. 1990). 

The difficulty in all such attempts originates in the relation 
between energy estimation with air showers and the mass of 
the primary nucleus. Generally showers initiated by heavy nu- 
clei develop and are absorbed faster in the atmosphere. Heavy 
nuclei thus produce showers of smaller size at the observa- 
tion levels of all existing experiments. The effect is stronger for 
inclined showers, which have to penetrate through larger atmo- 
spheric thickness. In this paper we model the shower size that 
the primary cosmic ray flux generates in the Akeno detector 
(Nagano et al. 1984), using both vertical showers and slanted 
showers; having performed such a fit, we then reconstruct the 
all particle spectrum. 

The paper is organized as follows: First we describe our in- 
put data and the most important parameters; then we present 
the data we wish to fit (Akeno) and the shower code used 
to calculate the observable parameter; our successful fits are 
shown with a discussion of possible errors in another section; 
we conclude with a discussion of the newly derived all particle 
spectrum and an outlook at possible next steps of interest. 



2. Input parameters 

Various cosmic ray experiments have given data about the 
chemical composition near TeV energies and somewhat be- 
yond, up to 100 TeV, most notably from the Chicago Group 
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(Grunsfeld et al. 1988, Miiller 1989, Swordy et al. 1990, and 
Miillcr et al. 1991), the JACEE experiment (Parnell et al. 1989, 
Burnett et al. 1990, and Asakimori et al. 1991a, b ), Simon 
et al. (1980) and Engelmann et al. (1990). As an illustrative 
example we show in Fig. 2 the Oxygen data, where we have 
averaged existing data from all four experiments when inde- 
pendent measurements exist in the same energy range. The 
data for the different experiments are quite consistent with 
each other. Energy ranges covered only by a single experiment 
are not represented. 



wind shock accelerated particles. We will therefore explore a 
range for the spectral index of the wind particles below the 
knee close to 2.67. 

2) The spectral difference at the knee. Again, simple theory 
(CR I) gives for this number the value 0.30, see above. The 
data on the spectrum above the knee suggest a spectral index 
of 3.0 ± 0.1. We will try a small range for the spectral index of 
the particles beyond the knee close to 3.00. 

3) The particle energy at the cutoff for the Sedov phase 
explosion acceleration. The Jacee data suggest for this num- 
ber approximately 100 TeV for Hydrogen. Our model firmly 
requires that this energy scale with the charge Z of the nucleus. 

4) The particle energy of the knee. Observations suggest 
an energy overall of about 5 10 3 TeV. Our theory requires that 
this energy also scale with nuclear charge Z, and if the knee is 
indeed dominated by heavy nuclei this energy for Hydrogen is 
likely to be above 200 TeV. 

The division of the abundances between the Sedov phase 
accelerated particles and the wind explosion accelerated parti- 
cles can be made using existing data, especially from JACEE. 



3. The data to be fitted and the shower development 
Monte Carlo code 



Fig. 2. The Oxygen spectrum derived from all existing data by bin- 
ning and then averaging with the weight of the error bars. A fit is 
shown to the resulting average spectrum showing a —2.64 spectral 
index. 



Data clearly demonstrate a power law behaviour with a 
somewhat flatter spectrum than the low energy index of about 
2.75. Oxygen has a spectrum of about 2.64, corresponding 
well to the argument made in CR I, that supernova explosions 
into existing stellar winds produce not only higher particle 
energies, but also spectra flatter than the lower energy particles 
from Sedov phase explosions into the interstellar medium. We 
use all these known chemical abundances as a given input, 
combining the elements into six groups, and the sources into 
three different sites as discussed above. The six element groups 
which we use are: a) Hydrogen, b) Helium, c) Carbon, Nitrogen 
and Oxygen, d) Neon to Sulfur, e) Chlorine to Manganese, 
and f) Iron. We assume that the galactic cosmic rays cut off 
at 10 Z GeV and that the extragalactic component is nearly 
all protons. In paper CR II (Biermann and Cassinelli 1993) 
we argue that this is in fact suggested by a number of stellar 
observations. 

The main uncertain parameters are then the following: 
1) The spectral difference between the Sedov phase explo- 
sion particles and the wind explosion particles. Simple first 
order theory gives 0.08 for this parameter, the difference be- 
tween 2.75 and 2.67, as described above. Arguments given 
in paper CR I and in Silberberg et al. (1990) suggest that 
the heavier elements in the TeV range come dominantly from 
shock acceleration in winds; therefore, we ought to look for 
a slight difference in spectrum between Hydrogen and heavy 
nuclei. For example, for Oxygen we find 2.62 ± 0.03 and for 
Hydrogen 2.74 ± 0.02, which corresponds to a spectral differ- 
ence 0.12 ± 0.04, consistent with our theoretical expectation 
for the difference between the Sedov phase particles and the 



The Akeno air shower experiment measures the size of showers 
arriving at the detector at different zenith angles. The result- 
ing shower size spectra are converted into constant intensity 
curves, i.e. groups of showers having the same experimental 
rate at different angles, thus initiated by primary nuclei of the 
same energy. 



Fig. 3. Average shower profiles for a primary particle energy of 10 
GeV for Hydrogen, Nitrogen and Iron. The two different slant depths 
which wc use arc indicated. 
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The angular dependence of the constant intensity curves 
is then used to extrapolate to size at shower maximum, from 
which an estimate of the primary energy that corresponds to 
a given intensity can be made. This fairly model independent 
procedure is still influenced by the chemical composition of 
the cosmic ray flux, especially if it has a strong energy depen- 
dence. The reason for this is illustrated in Fig. 3 which shows 
the development of showers with the same primary energy but 
different nuclei. We have decided to use directly the shower size 
spectra to reconstruct the cosmic ray spectrum under the as- 
sumption of different chemical compositions. We do this simul- 
taneously for two zenith angles, vertical (secanO — 1.0) and 
for secanO = 1.2, i.e. one moderate slant angle. At these angles 
the statistical errors are small. The shower size is expressed 
by N e - the total number of charged particles at observation 
level. The flux of showers with iV e = 10 6 at secanO = 1.2 is 
smaller than the vertical one by a factor of 7. The shapes of 
the size spectra are also different, as seen in Fig. 4. 

We use a simple parametrization (Gaisser, 1979) to cal- 
culate the size at the two slant depths (920 g/cm 2 for 
secanO = 1.0 and 1104 g/cm 2 for secanO = 1.2). While 
the original formula is for a constant inelastic cross-section, 
in this estimate we use a proton inelastic cross-section on air 
growing with the energy as (log E p ) 1& (Gaisser et al., 1987). 
For all nuclei heavier than Hydrogen we assume superposition, 
i.e. that all constituent nucleons interact independently in the 
atmosphere. This assumption is known to represent correctly 
the average shower size and underestimate its fluctuations. 

Our technique is intermediate between numerical integra- 
tion and Monte Carlo. We step through the energy spectrum of 
each chemical component in small (10 01 ) steps and calculate 
the shower size generated at the two depths by a nucleus of 
energy E. To estimate the fluctuations in N e , which are very 
important in view of the steep energy spectra, we calculate 100 
showers at each step, sampling the interaction depth for each 
nucleon. The resulting N e are binned with the appropriate 
weights. 

Although generally correct, the procedure used is not exact. 
A proper comparison of the cosmic ray spectrum and compo- 
sition model with experimental data should account fully for 
the fluctuations in the shower development (i.e. use a non- 
superposition nuclear fragmentation model as in Engel et al., 
1992) and those induced by the detection technique. This is 
a long term project which requires a close collaboration with 
the experimental groups. The technique used here, which is 
very sensitive to the cosmic ray spectrum and composition 
as we show further down, is fully sufficient as a first pass in 
the analysis of air shower data in terms of rapidly changing 
composition. 

The most important point here is to remember that shower 
size is reduced for particles at a given energy with higher 
nuclear mass: e.g. an Iron particle and a proton of the same 
energy will produce very different shower sizes, with the shower 
size of the Iron nucleus being much smaller (see Fig. 3) . 



4. The model fit 

We present the fits of the shower size distributions in Fig. 4, 
where the bands represent the experimental errors. 



Fig. 4. The fit of our model to the Akeno shower size data, both for 
vertical showers as well as for a slanted shower. Shown are the data 
as small boxes with the shading indicating the 1 sigma error range, 
as well as the model as a bold histogram. 

The comparison demonstrates that the fits are well within 
the acceptable error ranges. The model has the following pa- 
rameters: 

The cutoff energy (exponential cutoff) for the Sedov phase 
accelerated particles is 120 TeV for protons; the spectral index 
of the wind component is 2.66 below the knee and 3.07 above 
the knee; this steepening occurs at 700 TeV for protons and 
is rigidity dependent. The ratio of the polar cap component 
to the steeper wind component is 1 for heavy nuclei at the 
bend. Although the polar cap component does not contribute 
to shower sizes significantly below or above the knee, it is 
essential for reproducing the sharp break at iV e = 10 6 . The 
reproduction of this shower break, as well as the observed 
continuous flattening of the spectra of all heavy nuclei requires 
the introduction of a very flat (7 = 2) acceleration component, 
i.e. the polar cap component. We see that all four major 
parameters are close to their expected values. In Fig. 5 we 
show the contributions of the six chemical composition groups 
which we have distinguished, H, He, CNO, Ne-S, Cl-Mn, and 
Fe. We emphasize that all three critical particle energies, the 
cutoff energy for the Sedov phase accelerated particles, and 
both the bend energy and the cutoff energy for the wind 
accelerated particles are proportional to the nuclear charge Z, 
and are thus in our model fit 120 Z TeV, 700 Z TeV, and 10 5 Z 
TeV, respectively. 
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We also note that the code we use for reproducing the 
shower size distributions itself depends rather strongly on 
model fits to accelerator data to describe high energy in- 
teractions of nuclei; this is especially important for slanted 
showers where the tail of the shower profile is measured. Thus, 
the particle physics used introduces an additional uncertainty 
which is difficult to assess. 



5. Conclusions 



In Fig. 6 we show the all particle spectrum which results from 
our analysis. We compare it with the conventionally derived 
all particle spectrum, which does not explicity account for 
changes in composition. 



Fig. 5. The various contributions from the six element groups we 
distinguish: H (solid), He (dots), CNO (dash), Nc-S (long dash), 
Cl-Mn (dash dot), and Fe (long dash dot). 



ft is clear that within the framework of our picture the 
knee feature in the overall spectrum can only be due to the 
addition of the polar cap component. Without it the generated 
size spectrum cannot have a sharp bend in any model similar 
to ours, with rigidity dependent features and composed of 
different nuclear components. 

In presenting the fits we could have achieved smoother 
histograms by applying an exponential cutoff also to the polar 
cap component. We prefer to present the histograms in their 
present form, because a closer inspection now allows the reader 
to identify the component that causes the bend in the shower 
size spectra. In reality every single nuclear component has its 
own rigidity dependent behaviour which of course generates 
both smoother particle spectra and smoother size spectra. 

We note that the shower size distribution is very sensitive to 
changes in the parameters. Relatively small changes (by more 
than 20%) of the cutoff and bend energy, as well as in the 
composition of the wind component, affect the calculated size 
spectra so strongly that they become inconsistent with data. 
Together with the direct measurements in the 1 TeV region 
the shower size distributions constrain the models within a 
small parameter space. Fitting the shower size distribution at 
different zenith angles is critical despite its fairly large error 
bars. 



Fig. 6. Comparison of our all particle spectrum with the conven- 
tionally derived spectrum, as well as the various contributions from 
the different element groups. Here the symbols denote the following 
experiments: open circles, Akeno; open squares, Haverah Park; open 
triangles, pointing down, Yakutsk; open triangles, pointing up, Tien 
Shan; open hexagons, Fly's Eye; full squares, Proton 4; full circles, 
Jacee. The various references are given in detail in Stancv (1992) 
and in Hillas (1984). 

Several properties stand out: 

1) The correct spectrum is lower in the knee region and 
beyond than the conventionally derived spectrum; the factor is 
about 2. The knee itself is not as sharp but the change of the 
spectral index takes place at approximately the same energy. 

2) The composition changes rapidly in the knee region, 
becomes increasingly heavier, and in our current interpretation 
is dominated by the Neon group above the knee. One should 
note that the derivation of the exact charge of the dominating 
group depends strongly on the cascade model and that many 
observable parameters are not very different for Ne- and Fe- 
induced showers. The data will, however, be inconsistent with 
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a component lighter than Ne dominating the cosmic ray flux 
above 10 7 GeV. 

3) The proton flux shown in Fig. 6 at energy above 10 8 Gev 
represents the extragalactic component. Because the wind com- 
ponent has a maximum energy of 10 8 GeV for protons (2.6 10 9 
GeV for Fe) the calculated shower size spectra would become 
too steep for sizes above 10 7 without this extra component. 
Although the exact shape of the extragalactic spectrum needs 
further adjustments, its magnitude is limited by the compar- 
ison with the size spectra. This implies an extragalactic flux 
(presumably protons) that is consistent with the independently 
derived estimate by Gaisser et al. (1993) using the Fly's Eye 
data; a full comparison with an extragalactic model (Rachen 
and Biermann, 1992, 1993: paper UHE CR I) has been done 
in Rachen et al. (1993, paper UHE CR II). 

4) The comparison with existing direct data on the chemical 
composition is shown in Fig. 7 for the particle energy range 
1 to 100 TeV. Hydrogen, Helium and Iron groups data are 
shown. 



Fig. 7. A comparison of the measured element abundances for Hy- 
drogen, full squares, Helium, full circles, and Iron, full stars, with our 
model, again together with the all particle spectrum. The sources of 
the data are given in the first paragraph of section 2. 

The comparison demonstrates that our curves fit all ex- 
isting data quite well, including the trend, observed in heavy 
nuclei, to exhibit a flattening of the spectrum at the approach 
to the knee. The seeming inconsistency with the highest energy 
direct data for iron is possibly due to the fact the the JACEE 
experiment (where these points come from) has presented its 
measurements for nuclei with Z > 17. It is difficult to judge 
what the absolute normalization of the cosmic ray flux be- 
yond the knee is. The presentation of the spectrum in the form 
used here (S 2 75 dN/dE) tends to exaggerate differences in the 
absolute normalization. A relatively small error of the energy 
determination in air showers (typical error of 20%) changes the 
normalization of the absolute flux by a large amount (65%). 



We conclude that the most stringent test yet of the model 
proposed in the earlier papers of this series (CR I, CR II, CR 
III, UHE CR I, UHE CR II), which give predicted spectral 
shapes for all three sites of origin for cosmic rays is successful. 
Several detailed conclusions can be made: 

1) The abundances of the cosmic rays from explosions into 
stellar winds do not require nor imply any admixture from the 
heavy elements newly produced in the star exploding. All the 
heavy elements already present in the stellar wind prior to the 
explosion participate in the feeding of the accelerated particle 
population, and no other additional source is required. 

2) The nuclei accelerated in the polar cap region are essen- 
tial for a successful fit of the shower size spectra. Although the 
total energy carried by such nuclei is not more than 1/100 of 
the total wind component, they contribute a major fraction in 
the knee region. 

3) There is no requirement for other cosmic ray sources, 
either from spectral arguments or from abundance arguments; 
thus pulsars and compact X-ray binary systems may accelerate 
lots of particles, but they need not play a dominant role out 
in the typical interstellar medium. 

4) The large fraction of heavy nuclei inferred from the Fly's 
Eye data (Gaisser et al. 1993) requires the acceleration of heavy 
element nuclei out to at least 3 10 9 GeV. This implies (see CR 
I and II) that the winds of the stars that explode as supernovae 
have magnetic fields at least as strong as 3 Gauss at a fiducial 
distance from the star of 10 14 cm. Since this limiting particle 
energy is already derived by using the spatial limit given by 
the condition that the Larmor radius of the particle fit into 
the space available, there is no other way than indeed having 
these high fields in the stellar winds of massive stars. Should 
a more sophisticated analysis of the Akeno and Fly's Eye 
data require a galactic Cosmic Ray contribution of iron nuclei 
even at 10 10 GeV (see the discussion in Rachen, Stanev and 
Biermann 1993, UHE CR II), then a correspondingly higher 
magnetic field strength is required, with a lower limit of near 
10 Gauss at 10 14 cm from the star. 

The next steps possible and maybe desirable with better 
data or a more complete analysis of the existing data by the 
experimenters are the following: 

1) One must refine the cascade code to describe more ac- 
curately the high energy range and the highly slanted showers. 
Several different valid particle physics extrapolations to ul- 
trahigh energy should be used for an estimate of the errors 
caused by the uncertainty of the particle physics input. In 
addition, a non-superposition fragmentation model should be 
used to account more correctly for the shower development 
fluctuations. 

2) Also, one might combine all such source based calcula- 
tions with a proper propagation model in the Galaxy; after all, 
many of the heavy nuclei break up due to spallation. The task 
then would be to fit also the well known chemical abundances 
of the cosmic rays at GeV energies. 
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